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FOREWORD 


This  report  was  prepared  In  the  Mechanics  and  Surface  Interactions 
Branch  (AFWAL/MLBM) ,  Nonmetal 11c  Materials  Division,  Materials  Laboratory, 

Air  Force  Wright  Aeronautical  Laboratories,  Wrlght-Patterson  Air  Force 
Base,  Ohio.  The  work  was  performed  under  the  support  of  Project  Number 
2307,  "Nonmetalllc  Structural  Materials",  Task  Nunber  2307P2,  "Composite 
Materials  and  Mechanics  Technology  •  " 

In  this  report,  an  automated  composite  laminate  sizing  technique 
Is  presented,  which  optimizes  for  minimum  weight.  The  technique  can  be 
coded  for  a  microcomputer  and  a  listing  Is  given  for  the  Tlmex-Slnclalr  1000. 
The  program  Is  Interactive  and  easy  to  use.  Ply  ratios  are  optimized  for 
point  stress  under  multiple  Independent  loads. 

This  program  Is  available  on  cassette  tape  and  can  be  obtained  by  sending 
a  blank  15  or  30-minute  tape  to  AFWAL/MLBM,  Wrlght-Patterson  AFB,  Ohio  45433 
and  referencing  this  report. 


TABLE  OF  CONTENTS 


SECTION  PAGE 

I  PROGRAM  DESCRIPTION  I 

II  GENERAL  INSTRUCTIONS  2 

III  METHOD  6 

APPENDIX  A:  Representative  Results  13 

APPENDIX  B:  Subroutine  Descriptions  19 

APPENDIX  C:  Variable  List  23 

APPENDIX  0:  Notes  on  Sinclair  Basic  29 

APPENDIX  E:  Listing  31 


v 


SECTION  I 


PROGRAM  DESCRIPTION 


CLASS  (£omposite^aminate^utomated  £iz1ng  for  Strength)  is  an 
interactive  optimization  program  designed  to  run  on  a  small  microcomputer. 
The  listing  presented  here  is  for  a  Sinclair  ZX81  or  a  Timex-Sinclair  1000 
microcomputer  with  a  16K  memory  expansion.  The  version  of  Basic  is  standard 
enough  that  translations  to  other  microcomputers  is  possible. 

The  program  will  find  a  minimum  thickness  laminate  which  will  not  fail 
under  any  of  the  load  conditions  entered.  Ply  orientations  are  chosen  by 
the  user.  The  program's  capability  in  handling  multiple,  independent, 
loads  could  be  useful  for  loads  which  change  with  time  or  for  situations 
where  there  is  uncertainty  in  calculating  the  loads.  As  the  program  is 
currently  dimensioned,  four  independent  load  combinations  and  18-ply 
orientations  can  be  entered. 

Only  point  stresses  are  considered,  thus  the  program  optimizes  the 
laminate  only  at  one  point  in  the  structure.  Furthermore,  the  proqram 
assumes  in-plane  loads  only  and  no  out-of-plane  deflections.  This  implies 
a  symmetric  laminate,  but  stacking  sequence  is  not  a  factor  in  the  program. 
The  layer  thicknesses  generated  by  the  program  are  the  total  and  must  be 
divided  by  2  to  get  the  halves  of  a  symmetric  laminate. 

No  knowledge  of  optimization  techniques  is  needed  to  run  the  program 
and  very  little  knowledge  of  laminate  plate  theory  is  needed.  In  addition, 
material  properties  for  five  common  advanced  composites  are  stored  in  the 
program,  or  the  computer  can  ask  for  new  properties  through  prompts. 
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SECTION  II 


GENERAL  INSTRUCTIONS 

The  Timex  1000  or  Sinclair  ZX81  manual  Includes  taoe  loading 
Instructions.  Because  CLASS  takes  so  long  to  read  (approximately  7 
minutes).  Its  a  good  Idea  to  test  tape  player  volume  level  with  a  short 
one  or  two  line  program  to  see  If  all  is  well.  Load  the  tape  using  "CLASS" 
as  the  name.  If  the  tape  loads  properly,  it  will  automatically  begin 
execution.  An  example  of  the  video  prompts  and  appropriate  responses  are 
given  in  this  report.  As  a  number  is  entered,  it  appears  at  the  bottom  of 
the  screen.  The  number  can  be  changed  using  the  delete  key  (shifted  zero). 
Once  ENTER  is  pressed,  there  is  no  way  to  change  entries  until  the  end  of 
the  input  sequence  when  the  nrogram  asks  for  corrections.  If  there  are 
mistakes,  answering  "Y"  to  this  query  will  restart  the  sequence.  A 
subroutine  is  included  for  making  a  hardcopy  of  input  and  results  if  a 
printer  is  available.  At  the  completion  of  the  routine,  after  all  results 
have  been  displayed,  the  program  will  restart  itself. 

Run  times  can  be  quite  long.  They  range  from  a  minute  for  a  2-layer 
laminate,  to  an  hour  for  an  18-layer  laminate  subject  to  multiple  loads. 

Much  of  the  information  included  in  this  report  is  intended  for  those 
who  wish  to  understand  and  modify  the  program  and  is  not  needed  to  run  it. 


Display 

1 )  T300/5208 

2)  BORON/5505 

3)  AS/3501 

4)  SCOTCHPLY  1002 

5)  KEVLAR  49/EPOXY 

6)  AVAILABLE 

7)  NEW 

SELECT  MATERIAL 

HOW  MANY  PLY  ANGLES 

ENTER  PLY  ORIENTATION  (DEGREES) 

PLY  1  = 

PLY  2  = 

ENTER  NUMBER  OF  INDEPENDENT  LOADING 
CONDITIONS 

ENGLISH  OR  SI  UNITS  (E/S) 

LOADING  CONDITION  1 
N1  = 

N2  = 

N6  = 

ENTER  ENGINEERING  CONSTANTS  IN  GPA 
EX  *  ? 

EY  *  ? 

VX  «  ? 

ES  =  ? 


Keyboard  Response  and  Cownents 

(Tapes  will  automatically  start  after 
loading) 

(6  Is  an  available  slot  for  a  user 
defined  material.  Entering  7  allows 
new  properties  to  be  placed  in  the  slot. 
These  properties  can  be  used  in  sub¬ 
sequent  runs  by  entering  6) 


7  ENTER 
2  ENTER 

0  ENTER 
90  ENTER 

1  ENTER 
E  ENTER 

4000  ENTER 
1E3  ENTER 
^  ENTER 

(either  scientific  or  explicit  notation 
may  be  used  to  enter  numbers) 

181  ENTER 
10.3  ENTER 
.28  ENTER 
7.17  ENTER 
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1500  ENTER 


ENTER  STRENGTHS  IN  MPA. 

X  (TENS.)  =  ? 

X  (COMP.)  =  ? 

Y  (TENS.)  =  ? 

Y  (COMP.)  =  ? 

SHEAR  =  ? 

ENTER  MATERIAL  NAME  (<  15  CHARC) 


CORRECTIONS  (Y/N) 


(blank  screen  for  100  seconds) 

TOTAL  LAMINATE  THICKNESS 
.029706293  IN 

2  ACTIVE  CONSTRAINTS  AFTER  2 
ITERATIONS 

PRESS  ANY  KEY  TO  CONTINUE 


ANGLE 

RATIO 

NO.  PLIES 

0 

0.7542 

4.55 

90 

0.2458 

1.48 

PRESS  ANY  KEY  TO  CONTINUE 
STRENGTH  RATIOS 

1  =  ULTIMATE  STRAIN:  >  1  *  SAFE 
PLY  LOAD! 


1500  ENTER 
80  ENTER 
246  ENTER 
68  ENTER 

T30O/EPOXY(2)  ENTER 

(the  numbers  just  entered  are  simply 
those  of  T300/5208  with  double  the 
normal  transverse  strength.  If 
materials  1-5  had  been  selected,  all 
of  the  material  property  prompts  would 
have  been  skipped  by  the  computer) 

N  ENTER 

(this  simply  gives  the  user  a  chance 
to  restart  the  input  routine  before 
the  calculations  start) 


(actually  any  key  except  BREAK  which 
will  stop  the  program.  If  BREAK  is 
used,  the  program  can  be  restarted 
with  CONT) 


B 


0  1 

90  1.03 
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Press  any  key  to  continue 
HARDCOPY  (Y/N)? 

BAR  GRAPH  (Y/N)? 


B 

(useful  If  printer  Is  attached) 
N  ENTER 


1)  T300/5208 

2)  BORON/5505 

3)  AS/3501 

4)  SCOTCHPLY  1002 

5)  KEVLAR  49/EPOXY 

6)  T300/EP0XY(2) 

7)  NEW 

SELECT  MATERIAL 


(useful  If  many  ply  angles  are 
used.  Shows  relative  thickness 
of  layers) 

N  ENTER 

(program  automatically  restarts  at 
beginning.  Note  that  user  defined 
material  Is  now  number  6  and  can  be 
used  without  re-entering  the 
properties) 
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SECTION  III 


METHOD 


The  goal  Is  to  minimize  the  total  thickness  of  a  comoosite  laminate 
subject  to  failure  constraints  under  static  loads.  Specifically 
L 

Z  hi,  =  min.  where  L  =  number  of  layers 

k=l  K 

subject  to  hk  >  0 

and  -  1  £  0  where  h^  is  the  total  thickness  of 

all  the  plies  at  the  kth  orientation  (which  will  be  referred  to  as  a  "la^er" 
in  this  report).  The  failure  criteria  is  a  first  ply  failure  based  on  the 
Tsai-Wu  tensor  criteria  in  strain  space.  The  G's  are  transformed  to  the 
laminate  axis  from  the  k'th  layer's  orientation.  The  strains  are  associated 
with  the  N'th  loading  combination.  This  distinction  is  made  since  more 
than  one  independent  loading  may  be  considered.  For  the  definition  of  the 
G's  in  terms  of  experimental  strength  data,  see  reference  1. 

Stacking  sequence  is  not  included  in  this  formulation,  and  the  laminate 
is  assumed  not  to  bend  or  warp.  Therefore,  strains  and  loads  are  related  by 
N  =  |A|  e 

The  optimization  method  applied  is  a  modification  of  the  method  of 
2 

feasible  directions.  The  method  can  be  demonstrated  graphically  with 
2-dimensions,  i.e.  two  layers.  In  Figure  1  the  two  equalities 


e .£ .  +  G 
i  J 


(0) 

i 


Vj 


♦  G<90> 


-1=0 
e1  -  1  =  0 
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have  been  plotted  as  functions  of  and  h^^  for  the  single  loading 


condition  shown.  Any  point  above  and  to  the  right  of  these  two  curves 


is  feasible,  that  is,  failure  will  not  occur.  Points  to  the  left  and 


below  the  curves  are  infeasible.  Because  our  objective  function  (the 


sum  of  the  layer  thicknesses)  is  linear,  the  optimum  point  will  lie  on  one 


of  these  curves  or  the  intersection  of  multiple  curves. 


The  program  starts  by  finding  an  initial  feasible  point  (A)  which 

nn  a  f  aam  f  ha  ni  n  an  tha  lino  h[^)rh^^ 


lies  on  a  constraint  curve  farthest  from  the  origin  on  the  line  hl 


distance  from  the  origin  is  calculated  using  a  strain  ratio  method.  Along 


any  vector  which  passes  through  the  origin 
Li+1  _  ,i  .  * 


=  hk  *  $/so 


where  S  is  a  scalar  distance  and 


s0  -  Ck5,  (^)2],/2 


along  this  vector,  strain  can  be  found  using 


-  ei  Sr 


where  e?  is  a  component  of  laminate  strain  evaluated  at  S0.  Substituting 


into  the  failure  criteria  we  have 


G^’cV-S2  G(.ek,e°S 
TJ  1  J  0  +  1  0 

s2  s 


1  =  0 


To  ensure  the  calculated  point  lies  slightly  in  the  feasible  region  despite 


any  numerical  error,  the  program  sets  this  function  equal  to  the  negative 


of  a  small  number  (E^)  rather  than  zero.  Solving  this  equation  for  positive  S 


we  have 


C  _  -B  +  /(B-4AC) 

O  “  - XTT - 


I 
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where 


A  =  1  -  E1 

3  (8k) 

B  =  E  -  G.  £°S 

i=l 


Ji  1  o 


r  -  v  v  p(®|<)  0  °c^ 
C  =  E  I  -  G.  .K  £ .e  .S. 

1-1  j=l  1J  1  J  ° 


If  S0  lies  in  the  feasible  region  we  solve  the  above  equation  for  each  layer 
and  each  load  combination  then  take  the  smallest  resulting  S  as  the  one  that 
defines  the  boundary  of  the  feasible  region. 

The  next  step  in  the  optimization  procedure  is  to  establish  a 
direction  vector  which  will  point  away  from  the  constraint  A  lies  on  and 
is  parallel  to  the  plane  defined  by  Eh^  =  constant.  In  Figure  1,  this 
direction  is  shown  as  Z.  Finding  Z  first  requires  calculation  of  the 
gradient  of  the  active  constraint  evaluated  at  A.  Let 


G(0  )e(N)e(N)  +  G(e  )  (N) 
ij  i  j  i  i 


-  1 


where  k  and  N  correspond  to  the  layer  and  load  combination  of  the  active 
constraint.  A  constraint  is  considered  active  if 


Ck,N  -  "E2 

where  is  a  small  number.  Note  that  more  than  one  constraint  may  be  active. 
The  gradient  is  then  given  by 


*  c< 

L=1  1J 


8HW1  e(N)  +  (N)  ^jN) 
9h(_  J  i  3h 


)+  G 


(0k) 


9e 


(N) 


ah, 


A 


/V 

where  hL  is  a  unit  vector.  To  find  the  partials  of  strain  we  start  with 
the  basic  equation 
N  =  |A|  e 
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0  =  5^1*1  '*a3S7£ 


=  _ ! a- 1 !  3 


5-hT ,fl| 


(81) 

o(9,) 

12 

„(e  i) 
Q13 

11 

(ei) 

n(ei) 

n(8i) 

21 

q22 

Q23 

(ei) 

n(ei) 

(6i) 

g33 

13 

Q23 

The  gradient  vector  is  normalized  to  unit  length.  If  more  than  one 
constraint  is  active,  the  normalized  gradients  are  summed  together  and  the 
sum  is  then  normalized  to  one.  The  negative  of  the  gradient  will  point  away 
from  the  constraint,  into  the  feasible  region.  This  vector  is  now  projected 

A 

onto  the  plane  defined  by  the  unit  normal  n,  where 

"  =  k  hi 

The  projection  can  be  made  with  a  double  cross  Droduct 

/s  /S 

Z  =  n  x  (-Vc  x  n) 

With  a  vector  identity,  this  can  be  rewritten  as 
Z  =  (Vc  •  n)n  -  Vc 

Finally,  Z  Is  also  normalized  to  unit  length. 

Along  Z,  another  constraint  will  eventually  be  reached  (point  B  in 
Figure  1).  The  point  Is  found  Iteratively  by  a  bisection  technique.  Since 
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the  bisection  method  is  very  time  consuming,  the  constraint  line  is  only 
found  within  a  relatively  large  error  band.  What  we  are  really  interested 
in  is  a  point  approximately  midway  between  A  and  B,  which  is  C  in  the  figure. 
From  point  C,  the  strain  ratio  technique  is  used  to  analytically  calculate  D. 
Starting  at  D,  the  entire  procedure  repeats.  The  program  terminates  when 
the  distance  AB  or  CD  is  small  (say  1/10  a  ply  thickness)  or  the  magnitude 

■V  A  ^ 

of  Z  before  normalization  is  very  small  (implying  n  and  Vc  are  almost  parallel). 

In  some  cases,  hk  >_  0  constraint  may  be  reached.  When  this  happens, 
that  orientation  is  completely  dropped  from  further  calculations.  Thus, 
the  constraints  associated  with  a  zero  thickness  layer  cannot  effect  the 
results.  Once  an  orientation  reaches  zero  thickness,  it  is  never  reinstated 
in  later  iterations. 

Figure  1  shows  a  case  where  the  program  reaches  the  intersection  of  two 

constraints.  However,  simultaneous  failure  should  not  be  considered  a  criteria 

for  optimization.  Figure  2  shows  a  case  where  only  one  layer  approaches 

failure.  The  constraint  line  for  the  +45°  layer  is  completely  in  the 

[45]  [-45] 

infeasible  region.  The  line  h  +  hL  =  const,  has  been  included  to 
show  that  point  D  is  the  minimum  thickness. 
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APPENDIX  A 


REPRESENTATIVE  RESULTS 


i 

i 


Comments 


»»••*»»»•***»  »•»**»*•-•*-********* 
MATERIAL  T300/52O8 
TOTAL  THICKNESS  .000666667  H 
AFTER  3  ITERATIONS 


Nl=  1000000  N/M 
N2=  0  N/M 
N6  =  O  N/M 


ANGLE 

0 

90 

4-5 

-45 


RATIO 

1 

O 

O 

e 


NO.  PLIES 
5.33 

e 

e 

0 


STRENGTH  RATIOS 
1 -FAILURE:  HxSflFE 
PLY  LORD 1 

0  1 


********  ************************ 

?S3£ni?fcrIS5!Sr!!S07S39«ei2  M 

nrTER  4  rreRHTioNs 

ni»  1000000  n/m 

NS*  10OO0O0  N/M 
N6=  1000000  N/M 

Nls  2000000  N/M 
N8*  1000000  N/M 
N6*  0  N/M 

NU  2000000  N/M 

ns*  -leeeeoo  n/m 

NS*  e  N/M 

Nl*  0  N/M 
NS*  0  N/M 
N6*  -1500000  N/M 


ANGLE 

0 

45 

-45 

30 


RATIO 
0 . 3072 
O . 8335 
0.8471 
0.8122 


NO.  PLIES 
18.77 
14.27 
15.1 
12.97 


STRENGTH  RATIOS  _ 

1 -ULTIMATE  STRAIN:  > 1-SAFE 
PLY  LOAD1  LOADS  LOADS  LORD4 

0  1.848  1.819  8.398  1.078 

48  8.138  1.497  1.883  1.013 

>46  1.013  1.474  1.881  §.487 

90  1.307  1.894  1  1.87 


Demonstrates  program's 
capability  to  eliminate 
unnecessary  layers. 


Multiple  load  capability 
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Comments 


**•**********»*********-********* 
MATERIAL  T300/5208 
TOTAL  THICKNESS  .0055896009  H 
AFTER  3  ITERATIONS 

Nl»  2000000  N/M 
N2«  1000000  N/M 
N6*=  0  N/M 

N 1 =  1750000  N/M 
N2*  1250000  N/M 
N6>  433012.7  N/M 


ANGLE 

0 

90 

45 

-45 


RATIO 
0 . 3546 
0. 1583 
0 . 3064 
0 . 1807 


NO.  PLIES 
15 . 86 
7.08 
13.7 
8.08 


Note  that  this  */4  laminate 
and  the  tt/6  laminate  on  the 
next  page  give  the  same  total 
thickness . 

Second  load  has  the  same  magnitude 
as  the  first  only  with  the 
principle  axis  rotated  30°. 


STRENGTH  RATIOS 
1 =ULT IMATE  STR A IN :  >l=SflFE 
PLY  LOAD1  LOAD2 

0  1.261  1.123 

90  1.001  1.118 

45  1.034  1.307 

— 4S  1.2  1 


0.4 
0.36 
0.32 
0.28 
0.24 
0.2 
0.16 
0.12 
.08 
.04 

-90  -60  -30 


Bar  Graph.  Vertical  scale 
is  the  Dly  ratio;  horizontal 
is  orientation  angle. 
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Comments 


«4**»***-»i-«’**r.  **•*«***#•*■*•#■*•*■#•*•*•*■*■ 

HflTERIRL  T30O/5S08  „  .  , 

total  thickness  . 005589704 1  h  Same  loads  as  previous  example. 

AFTER  2  ITERATIONS 


Ml*  2000000  N /M 
NS=  1000000  N/M 
N5=  0  N.'M 
Nl  =  1750000  N.M 
N2  =  1250000  N/M 
N6s  4.33012.7  NsM 


ANGLE 

-60 

-30 

0 

30 


60 

00 


RATIO 
.099 
0 . 1564. 
0 . 24-26 
0.2434 
0 . 1584 
0 . 1002 


NO.  PLIES 
4.43 
7 

10.85 

10.88 

7.08 

4.48 


STRENGTH  RATIOS 
lsULTIMATE  STRAIN:  >1*=SAFE 
PLY  LORD 1  LOAD2 

-60  1.125  1.015 

-30  1 . 275  1 . 006 

0  1 . 254  1 . 103 

30  1 . 097  1 . 257 

60  1  1 . 275 

90  1.01  1.127 


0.3 


0.27 
0.24 
8.21 
0.18 
8. 15 
0.  12 
.09 
.06 
.03 


-90 
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18  layer  laminate 


*■*»•*»*****•  *••»*••***•*#*•**•** 
iflTCRiRL  Tsee/ssee 
TOTR‘_  THICKNESS  8.11247529  M 
RFTER  10  ITERRTIONS 

Nl=  4-0000000  N^M 
N2=  18000000  N/H 

N6  =  0  N/M 


Nl=  12500000  N/'H 
N2  =  32500800  N/M 
N6  =  12990000  N/M 


ANGLE 

RRT  IO 

NO.  PCXES 

-00 

.  809 

8. 13 

-70 

8 

0 

-60 

0 

0 

-50 

0 

8 

—4-0 

0 

0 

-30 

.0537 

48.35 

-20 

.099 

89 . 1 

-10 

0 . 1176 

105 . 82 

0 

0. 1059 

95 . 3 

10 

.  074-2 

66 .81 

20 

.0414- 

37 . 22 

30 

.  0252 

22 . 68 

4-0 

.0342 

30.3 

SO 

.0633 

56 . 98 

68 

.0963 

86 . 66 

70 

0. 1138 

102.41 

80 

0. 1028 

92.52 

90 

.0634 

57 

Second  load  is  of  the  same 
magnitude  as  the  first  but 
principle  axis  is  rotated  60°. 
These  loads  have  been  exaggerated 
by  an  order  of  magnitude.  This 
has  the  effect  of  making  the 
program  run  for  more  iterations. 


STRENGTH  RATIOS 
1 =ULT IHRTE  STRAIN: 


PLY 

CORD1 

LOROS 

-80 

1.023 

1.389 

-30 

1.449 

1.O04 

-20 

1.574 

1 

-10 

1.652 

1.015 

0 

1 . 649 

1.05 

10 

1.566 

1.11 

20 

1.44 

1.201 

30 

1.311 

1.325 

40 

1.201 

1.481 

50 

1 . 117 

1.648 

60 

1.058 

1.775 

70 

1.021 

1.801 

30 

1.004 

1.711 

90 

1.005 

1.554 

0.2 

9.18 

0. 16 

0.  14 

9.12 

0.1 

>i=SflFE 


.ee 

.06 

.04. 

.08 


-00  - 


Peaks  at  -10°  and  70° 


16 


Comments 


4********ir*5f*  ******************** 

;iflTERiHL  r3ee/5aes 

TQTRL  THICKNESS  .0110374-72  M  ,  ,  , 

=.fter  2  iterations  Same  loads  as  18  layer  example. 

•il=  4-000000  H/M 
N8=  100O000  N/M 

N6  =  0  N / M 


1250O00  N/H 
N2 -  3250000  M/M 
NS-  1239000  N/H 

ANGLE  RATIO  NO.  PLIES 

-10  e.5215  4-6 . 05 

70  0.4785  4-2.25 


Selecting  "peak"  orientations 
results  in  a  laminate  with  about 
the  same  total  thickness  as  the 
18  layer  case. 


strength  ratios 
1 =ULT iMRTt  STRRIN:  >1=SAFE 
PLY  LORP1  LORD2 

-10  1.364  1.007 

70  1  1 . 5S9 


******************************** 

MfUTEniRL  1300/5308 

TOTAL  THICKNESS  .023791632  N 

AFTER  3  ITERATIONS 

Nls  4000000  N/M 
N2=  1000000  N/H 

N6s  0  N/M  „ 

Same  loads  as  last  two  cases. 

Nl=  1250000  N/H 
N2  =  325O000  N/H 
N6«  1299000  N/H 


ANGLE  RATIO  NO.  PLIES 

O  0.2911  53.08 

90  0.7089  129.25 


STRENGTH  RATIOS 
1 .ULTIMATE  STRAIN:  >1=SAF£ 
PLY  LOAD1  LOAD2 

0  2.224  1 

90  1.32  1.018 


With  only  two  orientations, 
angle  sensitivity  is  important. 
Here  a  small  change  to  [0/90] 
has  resulted  in  a  laminate  twice 
as  thick  as  for  [-10/70]. 
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APPENDIX  B 


SUBROUTINE  DESCRIPTIONS 


i 

i 


500-760  CONSTRAINT  TEST 

Test  each  possible  failure  constraint.  If  a  constraint  is  violated,  set 
6$  *  "FAIL"  and  return.  If  no  constraints  are  violated,  set  6?  =  "PASS", 
make  a  list  of  active  constraints,  and  set  NC  (no  constraints). 


1000-1270  GRADIENT 

Find  the  gradient  of  the  constraint  identified  by  ply  P  and  load  N. 
Normalize  the  gradient  to  unit  length. 

1500-1760  STRAINS 

Given  a  value  of  S  find  |A|  =  } A^(  +  |Aj  xs.  Invert  [ A [  and  calculate 
laminate  strains  for  each  independent  loading. 


2000-2170  FORM  Ah  AND  A, 

L  *  L  L 

Ah  ii  =  E  A7  ••  =  E  ^iZk 

ij  K  z,ij  TJ  * 

Note  that  at  a  distance  S  along  Z,  | Af  =  |A^|  +  |Az|xs 


2200-2290  FORM  Q 

Convert  C  array  to  3x3  matrix  for  ply  II. 


2400-2498  FORM  G 

Convert  B  array  to  3x3  matrix  for  ply  II.  Also  form  vector  s  which  contai 
linear  components  of  failure  parameters. 


2500-2800  NEW  H 

Two  step  procedure.  First  find  approximate  distance  along  vector  Z  to 
next  constraint.  At  a  point  half  that  distance,  use  the  strain  ratio 


routine  to  decrease  the  total  thickness  with  constant  ply  ratios  until  a 
constraint  Is  found. 

3000-3350  NEW  DIRECTION 

Find  a  feasible  direction  vector  which  is  parallel  to  the  iso-thickness 
plane  and  leads  away  from  the  active  constraints. 

3500-3650  STRAIN  RATIO 

Along  a  vector  pointing  to  origin,  find  a  scalar  distance  from  current 
location  in  h  space  to  the  nearest  constraint. 

4100-4400  INITIAL  FEASIBLE  POINT 

Using  cn  assumption  of  equal  thickness  plies,  find  the  first  point  where 
no  constraints  are  violated.  Initialize  A^,  Az,  h,  strains,  and  constraint 
list. 

4500-4780  TRANSFORMATIONS 

For  each  ply  orientation,  transform  Q  and  6.  Store  the  results  in  B  and  C. 

5000-5495  INPUT 

Prompt  user  for  material,  angles  and  loads. 

5500-5658  OUTPUT 
Video  display  of  results. 

5660-5820  STRENGTH  RATIO 

Used  by  OUTPUT  and  HARDCOPY  to  print  out  a  list  of  strain  ratios  for  each 
loading. 
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5870-5950  PLY  RATIO 


Used  by  OUTPUT  and  HARDCOPY  to  print  out  a  list  of  ply  ratios  and  number  of 
plies. 


6500-6670  HARDCOPY 

If  printer  attached,  makes  a  printout  of  results. 


7000-7200  NEW  MATERIAL 

Prompts  user  for  new  material  prooerties. 


7500-7740  HISTOGRAM 
Generates  a  bar  graph. 


8000-8760  INVARIANTS 

Given  engineering  constants  and  strengths,  form  invariants. 
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APPENDIX  C 


VARIABLE  LIST 


Arrays 
A  (3,3) 

B  (18,9) 

C  (18,6) 

D  (3,3) 

E  (4,3) 

G  (3,3) 

N  (4,3) 

P  (3,3) 

Q  0,3) 

Vectors 
H  (18) 

R  (3) 

S  (3) 

T  (18) 

U  (5) 

V  (7) 

W  (18) 

X  (18) 

Y  (3) 

Z  (18) 


w 


a  matrix  for  current  value  of  thickness  vector 

contains  transformed  G's  (strength  parameters  in  strain  soace) 
for  each  ply  in  the  seouence  G-^ ,  G22.  G12,  Ggg,  G^g,  Ggg,  G^ , 

G2,  G3.  First  subscript  is  ply  number,  second  is  G  element. 

contains  transformed  Q's  (modulus  components)  for  each  ply 
in  the  sequence  Q^,  Q22,  Q] 2 »  $66’  $16’  Q26'  First  subscript 
is  ply  number,  second  is  Q  element. 

Ic 

I  Zk  where  Z  is  the  direction  vector 

strains  corresponding  to  each  independent  loading.  First 
subscript  identifies  load,  second  is  strain  component. 

(el’  £2’  e6} 

strength  parameter  matrix  for  a  given  ply  orientation 

loads.  First  subscript  identifies  independent  loading,  second 
is  load  element  (N-|,  N2,  Ng) 

A  inverse 

modulus  matrix  for  a  given  ply  orientation 


thickness  for  each  ply 
intermediate  results 

linear  strength  parameter  components  for  a  given  ply  orientation 

angle  of  each  ply  in  radians 

modulus  invariants 

strength  invariants  (strain  space) 

normalized  gradient  of  a  constraint 

normalized  sum  of  gradients  from  all  active  constraints 

Intermediate  results 

direction  vector 
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Scalars 

E2  -  defines  minimum  move  along  direction  vector  before  terminating 
program 

E5  -  defines  an  active  constraint  (if  G-jje.e.  +  G.e.  -  1  >  -E5  then 
constraint  active)  J 

E6  -  small  factor  included  in  strain  ratio  routine  to  guarantee  the 

point  found  is  slightly  in  the  feasible  region,  despite  numerical 
error 

C2  -  cos  2e 

C4  -  cos  40 

S2  -  sin  20 

S4  -  sin  40 

S  -  final  scalar  distance  to  be  moved 

51  -  point  in  feasible  region  in  bisection  routine 

52  -  point  in  infeasible  region  in  bisection  routine 

SREF  -  distance  to  origin  used  by  strain  ratio  routine 

SMAX  -  distance  along  direction  vector  to  first  h  =0  constraint 

511  -  units  conversion  Ib/in  -*•  N/m 

512  -  units  conversion  m  -*  in 
NPLY  -  number  of  ply  orientations 

NL  -  number  of  independent  loadings 
NC  -  number  of  active  constraints 
ITER  -  iteration  counter 
IMAX  -  iteration  limit 
M  -  identifies  material 

EX,  EY,  VX,  ES  -  engineering  constants  (reused  as  strength  parameters  in 
stress  space) 

XT,  XC,  YT,  YC,  SS  -  strengths 

QXX,  QYY,  Qxy,  Qs  -  modulus 

GXX,  GYY,  GXY,  GSS,  GX,  GY  -  strength  parameters  in  strain  space 
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TPLY  -  thickness  of  an  Individual  ply 
I,  0,  K,  L  -  loop  counters 
P,  II  -  ply  orientation  pointers 
N  -  load  pointer 

A.  B,  C,  CON,  DET,  E,  NORM,  TEST  -  intermediate  calculations 
Z  -  /NPLY*  where  NPLY*  is  number  of  ply  orientations  for  which  h^  f  0 

Strings 

C?  (10,2)  -  list  of  active  constraints,  identified  by  ply  and  loading 
P?  (6,19)  -  material  engineering  constants 

Y$  (6,24)  -  material  strengths 

M3  (6,15)  -  material  names 

R3  (4,2)  -  engineering  constant  labels  "EX",  "EY",  "VX",  "G" 

S$  (5,8)  -  strength  labels  "X  (Tens.)",  "X  (Comp.)",  etc. 

F3  -  flag  to  halt  program  when  =  "Fail" 

GJJ  -  flag  returned  from  constraint  test  routine  =  "Fail"  if  a  constraint 
is  violated 

K$  -  "Press  any  key  to  continue" 

U3  -  units  label  for  load 

V3  -  units  label  for  thickness 

E3  -  when  ■  "E"  English  units  are  desired 

A3,H?  -  query  responses 

Constants  entered  from  keyboard 

If  the  program  is  loaded  from  tape,  it  will  be  ready  to  run.  If  the 
program  is  keyed  into  the  computer,  certain  constants  and  dimension 
statements  have  to  be  entered  before  running  the  program.  These  have  not 
been  defined  in  the  program  in  order  to  save  memory.  When  the  program  is 
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SAVEd  on  tape,  the  constants  and  dimensions  will  also  be  saved.  Once  the 
constants  have  been  entered,  do  not  use  the  RUN  command  as  this  will  erase 
all  the  data. 


Dimension  Statements 


Dim 

A 

(3,3) 

Dim 

H 

(NPLY) 

Dim 

B 

(NPLY, 9) 

Dim 

R 

(3) 

Dim 

C 

(NPLY, 6) 

Dim 

S 

(3) 

Dim 

D 

(3,3) 

Dim 

T 

(NPLY) 

Dim 

E 

(NL ,3) 

Dim 

U 

(5) 

Dim 

G 

(3,3 ) 

Dim 

V 

(7) 

Dim 

N 

(NL ,3) 

Dim 

W 

(NPLY) 

Dim 

P 

(3,3) 

Dim 

X 

(NPLY) 

Dim 

Q 

(3,3) 

Dim 

Y 

(3) 

Dim 

c$ 

(10,2) 

Dim 

Z 

(NPLY) 

Dim 

Pi 

:  (6,19) 

Dim 

V9 

(6,15) 

Dim 

YS 

(6,24) 

Dim 

R  i 

(4,2) 

Dim 

Si 

(5,8) 

Where  NPLY  Is  the  number  of  ply  orientations  allowed  and  NL  is  the  number 
of  Independent  loads.  NPLY  =  18  and  NL  =  4  will  use  all  available  memory 
In  the  ZX81 . 

CONSTANTS 

M?(l)*"T300/5208" 

P?0)*"181 .  ,10.3,0.28,717" 

Y*(l  -"1500, 1500, 40. 0,246. ,68.0" 

M$(2)»"B0R0N15505" 

P?(2)*"204. ,18.5,0.23,5.59" 

Y?{2)-"1260, 2500,61 .0,202. ,67.0“ 

M*(3)-"AS/350r 
P?(3)*"138. ,8.96,0.30,7. 10" 

Y£(3)*"1447, 1447, 51 .7,206. ,93.0" 
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MS(4)*"Scotchp1y/1002" 

P$(4) *"38. 6, 8. 27 ,0.26, 4. 14" 
Y$(4)»"1062,610. ,31 .0,118. ,72.0" 

M*  5)*"Kevlar  49/Epoxy" 
P$(5)*"76.0,5.50,0.34,2.30" 

Y$( 5 )=" 1400,235. ,12.0,53.0,34.0" 

Let  SI1  *  175.1567 

Let  SI 2  *  39.37008 

Let  R*(l)  =  "EX" 

Let  R?(2)  =  "EY" 

Let  R*(3)  =  "VX" 

Let  R*(4)  =  "S" 

Let  SS(1)  =  "X(TENS. )" 

Let  S?(2)  =  "X(C0MP. )" 

Let  S?(3)  =  "Y(TENS. )" 

Let  S?(4)  =  "Y(C0MP. )" 

Let  S*(5)  *  "SHEAR" 

Let  K?  *  "PRESS  ANY  KEY  TO  CONTINUE" 

Let  E2  *  IE-5 

Let  E5  =  IE-1 

Let  E6  =  IE-6 

Let  TPLY  =  1.25E-4 

Let  IMAX  =  10 


27 


APPENDIX  D 


NOTES  ON  SINCLAIR  BASIC 


The  version  of  BASIC  used  on  the  ZX81  should  be  easily  translatable 
to  other  machines.  There  are  some  nonstandard  features  however,  which  may 
require  explanation. 

SLOW,  FAST  -  The  ZX81  uses  these  commands  to  control  whether  video  display 
is  continuous  or  goes  blank  during  computations.  They  can  be  ignored  for 
other  machines. 

A?  (a  to  B)  -  TO  is  the  string  slicing  command  and  replaces  the  standard  LEFT?, 
MID?  and  RIGHT?.  Note  that  string  slicing  is  used  to  define  material  properties. 
This  is  used  since  the  ZX81  lacks  READ,  DATA  and  RESTORE. 

LPRINT  -  Sends  string  to  printer. 

COPY  -  Sends  entire  video  display  to  printer 

PAUSE  40000  -  An  indefinite  pause,  broken  by  pressing  any  key 

Displays  in  the  program  are  designed  for  a  screen  that  has  21  lines  with 
32  characters. 
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APPENDIX  E 


LISTING 


MAIN 


11©  LET  ITER=1 

13©  oosub  see© 

138  PRINT  “CORRECTIONS  *?  CY/'N)  " 
134-  INPUT  8* 

13©  IF  R*="Y"  THEN  GOTO  130 
1*©  IF  H=7  THEN  GOSUB  700© 

15©  GOSUB  B©©© 

34-8  GOSUB  4-500 
84-©  GOSUB  4- 10© 

85©  GOSUB  3000 

86©  IF  F$="FBIL"  THEN  GOTO  550© 
89©  GOSUB  aS00 

3SS  Sr-nsKs*,!!™  soto  ssbb 

311  GOTOTis0I"RX  THEN  SOTO  SSBB 


S©S  LET  G$="Pfl5S" 

51©  LET  NC  =0 

580  FOR  P=1  TO  NPLY 

S25  IF  HCP)=©  THEN  GOTO  7S0 

58©  LET  lisp 

587  GOSUB  84-00 

S3©  FOR  N=1  TO  NL 

55©  LET  CON  =  —  1 

560  FOR  K  =  1  TO  3 

570  FOR  U=1  TO  3 

580  LET  CON=CON+G (K,U> »E(N, J) *E 
<N,K) 

590  NEXT  U 

©00  LET  CON=CON+S<K)#EtN,K) 

610  NEXT  K 

630  IF  CON >0  THEN  LET  G»  =  “F«IL" 
©4-©  IF  CON  <  — E5  THEN  GOTO  7©0 
©5©  LET  NC=NC+1 
©60  LET  C© (NC/ 1) zCHRt  P 
670  LET  C«<NC,8> =CMR©  N 

740  NEXT*n"FRIL’'  THEN  RETURN 
750  NEXT  P 
760  RETURN 


130  -  Call  innut 
132-136  -  Allow  user  a  chance 
to  change  input 
140  -  If  new  material  desired, 
branch 

150  -  Calculate  invariants 
242  -  Perform  transformations 
246  -  Find  an  initial  feasible 
point 

250  -  New  direction 
260,300  -  Halt  conditions, 
branch  to  output 
290  -  New  thickness  vector 
320  -  Loop  till  halt  condition 


CONSTRAINT  TEST 

525  -  If  nly  thickness  zero, 
ignore  constraints 
associated  with  it 
526-527  -  Set  up  G  matrix  for 
nly  being  tested 
560-610  -  Solve  con  =  G-.e-e-* 

GiEr’  u  1  J 

640-670  -  If  con  is  close  to 
zero,  identify  con¬ 
straint  as  active. 

C2  form  a  list  of 
constraints  in  terms 
of  ply  and  load 


GRADIENT 


ieie 
i«ia 
101* 
1020 
less 
1026 
1027 
1030 
104-0 
1045 
1050 
1060 
1070 
1060 
1000 
1100 
1110 
1120 
1136 
1150 
1155 
1160 
1170 
E  (N ,  K 
1180 
1190 
1192 
1194 
1200 
1210 
1220 
1250 
1260 
1270 


LET  NORM  — 0 

LET  I I =P 

SOSUB  2400 

FOR  L=1  TO  NPLY 

IF  HU.)  ==0  THEN  GOTO  1200 

LET  II=L 

GOSUB  2200 

FOR  0=1  TO  3 

LET  R  (J)  =0 

FOR  K  =1  TO  3 

nIxT*^’  *E  (N, 

NEXT  U 


FOR  0=1  TO  3 
LET  Y  CO)  =0 
FOR  K  =1  TO  3 
LET  V  CO)  =Y  CO) 
NEXT  K 
NEXT  O 
LET  U  CL)  =0 


♦P  (U,K)  *R  CK) 


K) 


FOR  0=1  TO  3 
FOR  K=1  TO  3 

LET  U  CL)  =U  (L)  +G  (J,K)  *  CY  fj)  * 
)  +E  CN,0)  *Y  CK)  ) 

NEXT  K 

LET  U  CL)  =U  CL)  +5  CO)  *Y  CO) 

NEXT  O 

LET  NORM=NORM+W  CL)  #U (L) 

NEXT  L 

LET  NORM  =S0R  NORM 
FOR  L=1  TO  NPLY 
LET  U  CL)  =U  CL)  /NORM 
NEXT  L 
RETURN 


1012-1014  -  Form  G  matrix  for 
designated  ply 

1026-1027  -  For  each  ply,  form  Q 
matrix 


1030-1070  -  R  *  -  fh  |A|  e 
1080-1130  -  Y  =  [A-1 ]  R 


-  Y 

8hi. 


Re, 


1150-1200  -  V(CON)  =  [G-,  (e. 

'J  i  dn^ 

3e .  a 

+  Jh!  E3>  +  Gf 


1210-1260  -  Normalize  v(CON) 


1510  FOR  1=1  TO  3 
1520  FOR  «J  =  I  TO  3 

1530  LET  Ea,J)=R(I(J)+Da,J)  *S 
1540  NEXT  O 
1550  NEXT  I 

1570  LET  DET=ECl,l)*E(a,a) *EC3,3 
) +2»E (1,2) *E (2,3) *E (1,3)  -E (2,2)  » 
E  Cl ,3)  *E  Cl ,3)  -E  Cl, 1)  *E (2,3)  *E  12, 
3)  -E (3,3) *E (1,2) *E (1,2) 

1580  LET  P Cl, 1) = CE (2,2) *E (3,3) -E 

(2.3)  *E  (2,3) ) 7DET 

1590  LET  P(2,2) =CE<1,1) *EC3,3> -E 

(1.3)  *E (1,3) ) /PET 

1600  LET  P Cl, 2) = (E (1,3) *E(2,3) -E 

(1,2)  *E  (3,3)  )  7DET 

1610  LET  P  (3,3)  =  (E  (1, 1)  »E  (2,2)  -E 

(1.2)  *E  (1,2)  )  -/DET 

1620  LET  P(l,3) =(E(1,2) *E (2,3) -E 

(2.2)  *E (1,3) ) 7DET 

1630  LET  P (2,3)  =  (E(1,2>  ME  (1,3)  — E 
(1,1) *E(2,3) ) /DET 
1640  LET  P (2, 1) =P (1,2) 

P  (3,2)  =P  (2,3) 

P (3,1) =P (1,3) 

1=1  TO  NL 
0=1  TO  3 
E (I,U) =0 
K=1  TO  3 

ECX,v!)  =E(X,wf)  +P(J,K)  *N( 


1650  LET 
1660  LET 
1600  FOR 
1690  FOR 
1700  LET 
1710  FOR 
1720  LET 
X,  K) 

1730  NEXT  K 
1740  NEXT  U 
1750  NEXT  X 
1760  RETURN 


FIND  STRAINS 

1510-1550  -  | E | i s  temporarily  the 
A  matrix  at  a  point 
along  the  Z  vector,  a 
distance  S  from  the 
current  position  in  h 
space 

1570-1660  -  Invert  A 

1680-1750  -  Solve  t  =  J A-1 J  N 

for  each  independent 
loading 
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2020 

2030 

204-0 

2050 

2060 

2070 

2080 

2090 

2095 

3100 

2110 

2120 

I) 

2130 

I) 

214.0 

2150 

3160 

2170 


2210 

2220 

2230 

2250 

2260 

2270 

2280 

2290 

2300 

2410 

34-20 

3430 

2440 

2450 

2460 

2470 

2480 

2490 

2492 

2494 

2496 

2498 


FOR  1=1  TO  3 
FOR  0  =  1  TO  3 
LET  9(1,0) =0 
LET  D  (I,  O) =0 
NEXT  O 
NEXT  I 

FOR  1=1  TO  NPLY 
LET  11=1 
GOSUB  2200 
FOR  0=1  TO  3 
FOR  K  =1  TO  3 

LET  R (U,K) =fi (U,K> +0 (U,K> *H ( 


FORM  Ah  AND  Az 


2120  -  Ah  is  the  A  matrix  at  the 
current  position  in  h  space 


2130 


A  .  . 

z,iJ 


NPLY 

Z 

k=l 


k 

q  .  z 

ij 


k 


LET  D (U,K)  =D (U,K)  +0 (U,K)  *Z  ( 

NEXT  K 
NEXT  O 
NEXT  I 
RETURN 


LET  0 
LET  0 
LET  O 
LET  O 
LET  0 
LET  0 
LET  0 
LET  0 
RETURN 
LET  G ( 1 


(1, 
(1, 
<  1 . 
(2, 
(2 
(3 . 
(  3  j 
(  3 , 


1)  =C (II, 1) 

2)  =C (11,3) 

3)  =C (11,5) 

2)  =C  (11,2) 

3)  =C (11,6) 

1)  =C (11,5) 

2)  =C (11,6) 

3)  =C (11,4) 


LET 

LET 

LET 

LET 

LET 

LET 

LET 

LET 

LET 

LET 

LET 


RETURN 


(1, 

(1, 

(2, 

(2, 

(2, 

(3, 

(3, 

(3, 

(1) 

(2) 

(3) 


1)  =B  ( 

2)  =B  C 

3)  =B  ( 

1)  =B  ( 

2)  =B  C 

3)  =B  ( 

1)  =B  ( 

2)  =B  ( 

3)  =B  ( 
=B  (II 
=B  (II 
=B  (II 


11,1) 

11,3) 

11. 5) 

11.3) 

11,2) 

11.6) 

11.5) 

11.6) 

11.4) 
,7) 
,8) 
,9) 


FORM  q 


2210-2290  -  Convert  C  array  into 

3x3  Q  matrix  for  a  ply 
designated  by  II 


FORM  G 

2410-2490  -  Convert  B  array  into 
3x3  G  matrix  for  a  ply 
designated  by  II 

2492-2496  -  Place  linear  terms  of 
G  in  vector  S 


33 


NEW  POSITION 


2501  LET  5MHX=1E10 

2502  FOR  1=1  TO  NPLY 

2503  IF  Z(I)<>0  THEN  LET  S=-H(X) 
✓Z  (II 

2504-  IF  S>0  RND  5<SMflX  THEN  LET 

5MflX=S 

2505  NEXT  I 

2507  LET  F*  =  "" 

2508  IF  SHRX  > 10  THEN  LET  F*=“FRI 

I  •« 


2510  LET  52=5HRX 

2512  IF  F*="FRIL"  THEN  RETURN 
2514-  LET  SI  =0 

2520  IF  NC=0  THEN  GOTO  2680 
2530  LET  S=S2 
254-0  GOSUB  1500 
2550  60SUB  500 

2560  IF  G*="FRIL"  THEN  LET  S2=S 
2570  IF  Gt^'PRSS"  THEN  LET  SI=S 
2585  IF  S1=SHRX  THEN  GOTO  2625 
2590  LET  S= (Sl+52) 72 

2600  IF  (52-51) <E2  RND  51=0  THEN 
LET  F*="FRIL" 

2605  IF  F*=”FRIL"  THEN  GOTO  2760 
2610  IF  51/- (52-51)  <4-  THEN  GOTO  2 
54-0 


2620  LET  5=5172 

2625  LET  SREF=0 

2630  FOR  1=1  TO  NPLY 

2640  LET  H  (I)  =H  (I) +Z  (I)  *5 

2650  IF  H ( I) <E3  THEN  LET  H(I)=0 

2660  LET  SREF=SREF+H (II  *H (I) 

2670  NEXT  I 

2680  LET  5=0 

2690  GOSUB  2000 

2700  GOSUB  1500 

2710  LET  SREF=SQR  SREF 

2720  GOSUB  3500  _  _ 

2726  IF  SREF  — S  <E2  THEN  LET  F*  =  "F 


RIL” 


2730  FOR  1=1  TO  NPLY 

274-0  LET  H  (I)  =H  (I)  *S/SREF 

2750  NEXT  I 

2760  LET  5=0 

2770  GOSUB  2000 

2780  G05UB  1500 

2790  GOSUB  500 

2800  RETURN 


2501-2505  -  Find  distance  along  Z 

to  first  h-=0  constraint 
2540-2610  -  Bisection  method  to  find 

distance  to  next  constraint 
2585  -  If  no  constraints  are  violated 
at  S  =  SMAX  then  stop  search 
and  use  that  point 
2620  -  Take  a  point  halfway  in 
between  constraints 
2625-2670  -  Uodate  h  vector  at  that 
point  and  calculate 
distance  to  origin 
2690-2700  -  Update  strains  at  the 
mi d-point 

2720  -  Use  strain  ratio  routine  to 
find  nearest  constraint  along 
a  line  from  the  midpoint  to 
the  origin 

2730-2750  -  Update  h  vector  to  new 
noint  near  constraint 
2770  -  Update  A  matrix 
2780  -  Update  strains 
2790  -  List  of  active  constraints 
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NEW  DIRECTION 


3005  LET  NORM-1 
3Ci0&  LET  Z-0 
3020  FOR  1=1  TO  NPLY 
3030  LET  xa.!  =0 
3940  LET  Z-Z+SGN  H  (  I  :> 

3050  NEXT  I 

3052  LET  Z  =1 /SOR  Z 

3060  IF  NC-0  THEN  GOTO  3225 

3070  FOR  J-l  TO  NC 

3072  LET  P=CODE  C-*(I,i) 

3074  LET  N=CODE  C*CI,2) 

3080  GOSUB  1000 

3140  FOR  0-1  TO  NPLY 

3150  LET  X  (J)  =X  CJ)  -U  (U) 

3160  NEXT  U 

3170  NEXT  1 

3180  LET  NORM  =0 

3190  FOR  0=1  TO  NPLY 

3200  LET  NORM=NORM+X (J) *X(J) 

3210  NEXT  O 

3220  LET  NORM  =SOR  NORM 
3224-  LET  TEST=0 

3225  FOR  1=1  TO  NPLY 

3226  LET  X  Cl)  =X(I)  /NORM 

3230  LET  TEST -TEST +X (  I )  *Z  *SGN  H( 
X) 

324-0  NEXT  I 
3250  LET  NORM  =0 
3260  FOR  1=1  TO  NPLY 

3270  LET  Z(I) =X(I) -TESTiZ*SGN  HC 
I) 

3280  LET  NORM=NORM+ZfI) *ZCI) 

3200  NEXT  I 
3202  LET  F$="" 

3204  IF  NORM  <  IE— 6  THEN  LET  F*  =  **F 
RIL" 

3206  IF  F#=“FRIL"  THEN  RETURN 

3300  LET  NORM=SOR  NORM 

3310  FOR  1=1  TO  NPLY 

3320  LET  ZCX)=Z(I)/NORH 

3330  NEXT  X 

3340  GOSUB  2000 

3350  RETURN 


3020-3052 

3070-3080 

3140-3320 

3225-3240 

3260-3290 

3292-3296 

3310-3330 


Initialize  x.  Z  is 
the  square  root  of  the 
number  of  non-zero 
thickness  plies 
For  each  active  constraint 
calculate  the  gradient 
of  the  constraint.  p 
and  N  identify  the 
constraint  to  the  gradient 
routi ne 

Sum  all  gradients  into 
x  and  normalize  x 
Test  is  the  dot  product 
of  x  and  the  unit  normal 
(h)  to  the  ::hi  =  const, 
plane 

— F 

Z  is  now  a  vector 
parallel  tG  = 
constant  plane  and 
pointing  away  from  the 
active  constraints  _ 

If  the  magnitude  of  Z 
is  clpse  to  zero  then 
Z  | |  h  and  a  local 
minima  has  been  reached 
normal ize  Z 


3510  FOR  P  =  1  TO  NPLY* 

3520  IF  H<P)=0  THEN  GOTO  3540 

3522  LET  II=P 

3524  GOSUB  2400 

3530  FOR  N=1  TO  NL 

3540  LET  B=0 

3545  LET  C=0 

3550  FOR  1=1  TO  3 

3560  FOR  U=1  TO  3  _ 

3570  LET  C=C-SREF*SREF*Q  tI,J) *E  C 
N,  I)  *E  CN  ,  U) 

3580  NEXT  U 

3590  LET  B=B-SREF*S Cl) #E CN, I) 
3600  NEXT  I 

3610  LET  SURL=(-B+S0R  tB#B-4*C# C 
1-E61 )  ) 7 (2*  C1-E6) > 

3680  IF  SURL  >  S  THEN  LET  S*SURL 
3630  NEXT  N 
3640  NEXT  P 
3650  RETURN 


STRAIN  RATIO 

3522-3524  -  Form  G  matrix 
3522-3610  -  For  each  constraint  solve 

Gij£ic.i(SREF)2 

2 


+  - 


G^e • (SREF) 


1  =  E6 


for  S 


3620  -  Take  smallest  vale®  of  S 
(closest  constraint) 
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4110 
4120 
4125 
4130 
4140 
4150 
4160 
4170 
4130 
4185 
4190 
4  330 
4340 
4350 
4360 
4375 
4380 
4390 
4400 


4510 
4520 
4530 
4540 
4550 
4560 
U  <3> 
4570 
'J  (3) 
4580 
4590 
4600 
) 

4610 

) 

4650 
4660 
4670 
4680 
U  (3) 
4690 
U  <3) 
4700 
4710 
4720 
) 

4730 

) 

4770 

4780 


INITIAL  FEASIBLE  PT. 


LET  Z  =  l/SOR  NPLY 

FOR  1=1  TO  NPLY 

LET  Z  <  I)  =Z 

LET  H  < I)  =  Z 

NEXT  I 

GOSUB  2000 

LET  5=0 

LET  SREF=i 

GOSUB  1500 

LET  5=0 

GOSUB  3500 

FOR  1=1  TO  NPLY 

LET  H  ( I)  =H(I)  *5 

NEXT  I 

LET  5=0 

GOSUB  2000 

GOSUB  1500 

GOSUB  500 

RETURN 


4110-4150  -  Initialize  the  h  vector 
to  an  arbitrary  DOint 
alonq  the  line  h.j=h2 

=  hNPLY 

4180  -  Establish  the  strains  at 
that  point.  Point  is  far 
enough  out  that  no  constraints 
are  violated  for  reasonable 
structures 

4190  -  Use  strain  ratios  to  find 
nearest  constraint 
4330-4350  -  h  vector  updated  so 
that  position  in  h 
sr>ace  lies  on  constraint 
4375-4390  -  Update  A  matrix  and 
constraint  list 


FOR  1=1  TO  NPLY 
LET  C2=COS  C  2  *T  ( X )  ) 

LET  C4-COS  ( 4  *T  ( X)  ) 

LET  S2=SIN  ( 2  *T  ( I >  ) 

LET  S4  =S IN  ( 4  *T  ( I )  ) 

LET  B  (I  ,  1)  =U  (1)  +C2*U  (2)  +C4* 

LET  B (1,2) =U (1) -C2*U (2) +C4* 

LET  B ( I , 3) =U (4) -C4*U (3) 

LET  B (I ,4) =U C5J -C4*U (3) 

LET  B  «I , 5)  =S2/2*y (2 J  +S4*V (3 

LET  B(I,6) =S2/2*U(2) -S4*U(3 

LET  B (Xj7) =U (6) +C2*U (7) 

LET  B(I,8)=U(6)-C2*U17) 

LET  B(X,9) =S2*U(7) 

LET  C < X , 1) =U < 1) +C2*U (2) +C4* 

LET  C ( X , 2) =U ( 1) -C2*U (2> +C4 * 

LET  0(1,3) =U(4) -C4*U(3) 

LET  C (1,4) =U (5) -C4*U (3) 

LET  C (I,SJ =S2/2*U (2) +S4*U <3 

LET  C (1,6) =S272*U (2) -S4*U <3 

NEXT  I 
RETURN 


TRANSFORMATIONS 


4560-4670  -  Transform  failure 


parameters  in 

following 

order 

B(I,1)=GU 

B ( I  ,5)=G16 

B(I,2)=G22 

B(I,6)=G26 

B(I,3)=G12 

B(I,7)=G] 

Bd.4)=G66 

B(I,8)=G2 

B(I,9)=G3 

4680-4730  -  Transform  modulus  in 
following  order 
C(I,1)=Q11  C(I,4)=Q66 

C(I,2)=Q22  C(I,5)=Q16 


C(I,3)=Q12  C(I,6)=Q26 


Note,  transformations  for  all 
orientations  calculated  and  stored 
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INPUT 


500P 

CL  5 

500" 

FflSI 

5090 

FOR  1=1  TO 

6 

5  100 

PRINT  X..  '  5 

51  10 

MFXT  I 

5  115 

PRINT  ‘7) 

NEU" 

5120 

SLOW 

5160 
XRI  ' 

PRINT  RT  IS 

;,S, "SELECT  HflTEP 

3  176 

INPUT  H 

5190 

CL.  5 

S200 

PRINT  "HOU 

MPiNY  PLY  ANGLES  7 

52  ie< 

1NPUT  NPLY 

5216 

ClS 

522  is 

r-'RlN'.'  ‘ENTER  PLY  OHISNTflTIU 

5.230 

FOR  I 

=  1  TCi  NPLY 

S24S 

PRINT 

"PLy  ‘‘.I. 

"  22  " 

5250 

INPUT 

Til) 

5260 

PRINT 

RT  I  .  7  ;  T  <  1 5 

5270 

L  IT  T 

{  I)  =T  (  I)  iP 

I  /  ISO 

5290 

NEXT 

I 

5295 

PRUSE 

120 

5300 

CLS 

53  10 

PRINT 

•ENTER  NUMBER  OF  INDE 

PE. -4  OR  NT 

LORDING 

CONDITIONS" 

5320 

INPUT 

NL 

33  30 

LET  U‘V  =  "  M/M" 

5  3  4-  0 

LET  U 

$  =  "  M " 

54  22 

PRINT 

“ENGLISH 

OR  SI  UNITS? 

1E,'£ 

•5  ’• 

5424.^ 

INPUT 

El- 

54-26 

IN’ 

F* 

=“E"  THEN 

LET  U*="  LBS 

54.28 

IF  E* 

=  ‘‘E“  THEN 

LET  U*  =  "  IN¬ 

3430 

FOP  I 

=1  TO  NL 

5440 

CLS 

5450 

J 

PRINT 

“LORDING 

CONDITION 

5455 

FOR  J 

=  1  TO  3 

3456 

LET  L 

=  J 

5457 

IF  U  = 

3  THEN  LET 

L=6 

546C 

PRINT 

N";  L, 

5465 

INPUT 

N  t  I  .  J) 

5470 

PRINT 

RT  J,4,N(I,J) ,U$ 

5475 

IF  E* 

a-E"  THEN 

LET  N i I,U) =N 

(  r ,  U  ) 

tSIl 

5480 

NEXT 

U 

5482 

PRUSE 

120 

34-85 

NEXT 

I 

5486 

CLS 

5490 

FAST 

54.95  RETURN 


5090-5160  - 


5210  -  NPLY 


5230-5290  - 


5330-5428  - 


5450-5480  - 


List  available  materials 


=  number  of  layers 


Enter  orientations 
computer  requires  angles 
in  radians,  so  convert 
degrees  to  radians 


Establish  proper  units 
labels 


Enter  loads.  If  loads 
are  in  Ibs/in,  convert 
to  N/m 
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OUTPUT 


5S®5 

5506 

5510 

5520 

5530 

5535 

553- 9 

554- 0 
NESS 
S54-4- 
554-9 
5550 

r  s  ” 

5555 

T  IONS 

5550 

5570 

5510 

5615 

5620 

5622 

5624. 

5626 

5626 

5630 

5632 

5634. 

5635 

5640 

5645 

5655 

5656 

5657 
:.65S 


LET  TEST  *0 
LET  E=1 

FOR  X»1  TO  NPLY  _ 
let  test  =TEST -fH  (I) 

NEXT  I  .  ___ 

IF  E*»"E“  THEN  LET  E*SI2 

PRINT  '"TOTAL  LAMINATE  THICK 

PRINT  TAB  10,TEST*E.i  U* 

PRINT 

PRINT  NC,"  ACT IUE  CONSTRAIN 

PRINT  "AFTER  " , ITER /  ITERR 

PRINT  RT  21j0i K ♦ 

PAUSE  40000 

CL  5 

(30  SUB  5870 
PRINT  RT  21,0..K* 

PAUSE  40000 
CLS 

GOSUB  5660 
PRINT  RT  21)0; K ♦ 

PAUSE  40000 
CLS 

PRINT  "HARDCOPY  (Y/N>7" 
INPUT  H* 

IF  H$= "Y“  THEN  GOSUB  6500 

PRINT  "HISTOGRAM  <Y /N) 7” 

INPUT  0*  _ _ 

IF  R$s"V"‘  THEN  GOSUB  7500 
GOTO  1 


5510-5535  -  Sum  for  total  thickness 
and  establish  units 
conversion 


5615  -  Call  ply  ratio  printout 
5626  -  Call  strength  ratio  printout 


56BG®PRINT  "STRENGTH  RATIOS  "" 
5665  PRINT  *"  1  =ULT IHRTE  STRAIN  :  > 

1=SAFE" 

5670  PRINT  "PLY- 
5672  FOR  Is=l  TO  NL 

5674  PRINT  RT  2  ,  I  *7  - 1;  " LORD "  I 

5676  NEXT  I 

5680  FOR  P  = 1  TO  NPLY 

5685  IF  M(P)=0  THEN  GOTO  5810 

366©  LETT  II=P 

5690  GOSUB  2400 

3700  FOR  N=1  TO  NL 

5702  LET  R  3t0 

5704  LET  B=0 

5710  FOR  U=1  TO  3 

5720  FOR  K  =  1  TO  3 

3730  LET  R=R+G  (U..K)  »E  (N,UJ  *E  (N,  K 
3740  NEXT  K 

5750  LET  8=B+S(J)*E(N,U) 

5760  NEXT  U 

5770  LET  R=(-B+SQR  (BiB+4iA) ) ✓ (2 
*R» 

5780  LET  RsINT  ( (R+ . S^1E3) *1E3) ✓ 
1E3 

5785  PRINT  RT  P+2,  N*7-l,fi 
5790  NEXT  N 

5800  PRINT  RT  P+3,0.;  U  <P) 

5810  NEXT  P 
5820  RETURN 


STRENGTH  RATIOS 


5700-5780  -  Solve  for  R  in 

r2  *  [vfj  r 

for  each  loading 


0 


PLY  RATIO 


3880  PRINT  "flNGLe";TOB  11; "RATIO 
".TAB  21, "NO.  PLIES’ 

3896  FOR  1*1  TO  NPLY 
5900  LET  U  < I )  * INT  (tT(I> 

.05) *10) /10 
5910  LET  p  a INT 
4.)  *1E4)  ✓  1E4 
3920  LET  B  =  INT 
*100) /100 
5930  PRINT  U(I) 

3 

394.0  NEXT  I 
3950  RETURN 


*180/PI+ 
(  (H  Cl)  /TEST* . 5/ IE 
C  (H  C I)  /TPLY+ . 005) 

;  TAB  11;  ft;  T AB  24.; 


Tor  each  orientation,  find  ply 
NPLY 

ratio  =  (h -/  T.  hj  and  number  of 
1  i=l  1 

plies  =  h.j/(ply  thickness) 


3510  LPRINT  '******************* 

*************"  .  „ 

3520  LPRINT  “MATSRIAL  ";M*(H) 
3525  LET  A$=5TR$  (TEST  *E)  +U* 

3530  LPRINT  "TOTAL  THICKNESS  ;  R 


5540  LPRINT  "RFTER  " 
RTIONS" 

5560  FOR  1*1  TO  NL 

5570  LPRINT 

3572  LET  L*1 

5574  IF  E**"E”  THEN 

5580  FOR  U=1  TO  3 

5582  LET  fi=«J 

5586  IF  R=3  THEN  LET 

559©  LPRINT  ''N'* ;  R;  "  * 

J$ 

5600  NEXT  U 
3610  NEXT  I 
3620  LPRINT 
3625  CLS 
3630  GOSUB  5870 
6640  COPY 
6645  CLS 
6650  GOSUB  5660 
6660  COPY 
6670  RETURN 


;  ITER;  "  ITER 


LET  L=1/5I1 


R*=6 

•• ;  N  ( I ,  vJ )  *L ; 


HARDCOPY 

5560-5610  -  Print  loads  with 
given  units 


5630  -  Call  ply  ratio 
5650  -  Call  strength  ratio 


NEW  MATERIAL 


7005  SLOU 

7006  CLS 

7007  LET  H=6 

7010  PRINT  "ENTER  ENGINEERING  CO 
NSTRNTS" 

7015  PRINT  "IN  GPR • 

7020  FOR  1*1  TO  4 

7030  PRINT  RT  1+1,0; R$ (1) ; "  *  7" 

7040  INPUT  R 

7050  PRINT  RT  1+1,4; fi 

7060  LET  P* (6, 1*5-4  TO  1*5-1) «ST 

R$  R 

7070  NEXT  I 
7080  CLS 

7090  PRINT  "ENTER  STRENGTHS  IN  M 

PR  " 

7100  FOR  1-1  TO  5 

7110  PRINT  RT  X+2, 0; S( (I) ; "b 

7120  INPUT  R 

7130  PRINT  RT  1+2, 10; R 

7150  LET  Y* (6, 1*5-4  TO  1*5-1) »ST 

R9  R 

>160  NEXT  I 
7170  CLS 

7180  PRINT  "ENTER  MRTERXRL  NRHE 
f  < 18  CHRRC . ) " 

7190  INPUT  M0tS) 

7198  FRST 
7200  RETURN 


7030  -  R£  contains  prompts  Ex, 

Ey»  Vx,  G 

7060  -  String  concatenation  to 
store  properties  in 
array 

7110  -  S3  contains  prompts  X(TENS.), 
X(C0MP. ) ,  Y (TENS . ) ,  Y(C0MP.), 
SHEAR 

7150  -  String  concatenation  to  store 
strengths  in  Y?  array 
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HISTOGRAM 


750B1LET  Z  =0 
75®2  CL5 

7505  FOR  1=1  TO  NPLY 

7510  IF  H (I)  /TE ST  >  Z  THEN  LET  Z=H 

(I) /TEST 

7520  NEXT  I 

7522  LET  DELTR  =  INT  (Z/9»100fl)  /I 
00 

7525  FOR  1=1  TO  NPLY 

7S3®  LET  Z<I)=INT  CH  (I)  *4- .  2/TEST 

/OELTR+.5J 

7540  NEXT  I 

7550  FOR  1=1  TO  NPLY 

7560  FOR  0=3  TO  Z(X)+2 

7570  FOR  K  =  —  1  TO  1 

7560  PLOT  INT  (U(I) *3/10+35.5) +K 

..  O 

7590  NEXT  K 
7600  NEXT  O 
7610  NEXT  I 
7620  FOR  1=1  TO  7 

7625  LET  U  =  INT  (I *4. 5 +.5) -2 

7626  IF  0=16  THEN  LET  0=0+1 
7630  PRINT  RT  21 , O ; -120+1*30 
7640  NEXT  I 

7660  LET  Z=11*DELTR 
7670  FOR  I=®  TO  18  STEP  2 
7680  LET  Z=Z -DELTR 

7690  LET  R=INT  ( (Z+.005) *1E2J /IE 

2 

7700  PRINT  RT  1,0; R 
7710  NEXT  I 

7720  IF  H*<>"Y"  THEN  PAUSE  40000 
7730  IF  H*="Y"  THEN  COPY 
7740  RETURN 


7505-7520  -  Find  largest  h^ 

7522  -  Vertical  scale  increment 
7530  -  Bar  height  in  terms  of  pixels 

7580  -  Plot  bars,  3  pixels  wide 

7620-7640  -  Horizontal  scale 

7660-7710  -  Vertical  scale 
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INVARIANTS 


3160  LET  EX=URL  P$(M,1  TO  4-)  *1E© 
8170  LET  EY»URL  P*(H,6  TO  9) *lEO 
8180  LET  UX«URL  P*(H,11  TO  14.) 
3190  LET  E9*URL  P$(M,16  TO  > *1E9 
3200  LET  XT *URL  Y*tM,l  TO  4>*1E6 
8210  LET  XCsURL  Y$rM,6  TO  SO  *1E6 
3220  LET  YT.URL  Y*CM,il  TO  14)  *1 

8230  LET  YC»URL  Y*<M,16  TO  19) *1 
£6 

3240  LET  5S«=URL  Y*(M,21  TO  )  *1E6 
3505  LET  N* Iy ( 1 -UX *OX*EY /EX ) 

3510  LET  OXX=N*EX 
3520  LET  OYYxN*EY 
3530  LET  QXY=N*OX*EY 
3535  LET  ©S*ES 

•3540  LET  U  (1)  x  (3*OXX+3*OYY+2*OXY 
+4*ES) sO 

3550  LET  U (2) = (OXX-OYY) /2 

3560  LET  U (3) s (QXX+QYY-2*QXY-4*E 

S)  ✓S 

8570  LET  U (4) * COXX+OYY+6*OXY-4*E 
S)  /'S 

8580  LET  U  (5)  =  (QXX+OYY -2  *OXY +4  *E 
5) 

3590  LET  EX  *1/1  XT  fXC ) 

3600  LET  EY=1/(YT*YC) 

3610  LET  ES  =  1/  (SStSS) 

3612  LET  FX=1/XT-1/XC 

3614  LET  FYsl/YT-l/YC 

3620  LET  EXY=-SQR  (EX*EY) ^2 

3630  LET  GXX=EX*OXX*9XX+2*EXY*©X 

X  *OXY  +EY  *QXY  *©XY 

8640  LET  GYY=EX*QXY*QXY+2*EXY*OX 
Y  *OYY  +EY  *©YY  fOYY 

3650  LET  GXY  =EX  *OXX iOXY +EXY  * ( OXX 

tOYY  +OXY  J-OXY  )  +EY  *OXY  *OYY 

3660  LET  GSS=ES*OS*OS 

3670  LET  GX=EX*OXX+FY*QXY 

3680  LET  GY=FX*OXY+FYtOYY 

3690  LET  U ( 1) = (3*GXX+3*GYY+2*GXY 

+4*G5S) /& 

3700  LET  U(a)«(GXX-GYY)/S 

3710  LET  U (3) = (GXX+GYY-2*GXY-4*G 

35)  /& 

3720  LET  0(4)* (GXX+GYY+6*GXY-4*G 

36)  /8 

3730  LET  O (5) * (GXX +GYY -2  *GXY +4  *G 
SS)  y-S 

3740  LET  U  (6)  =  (GX+GY)  /-2 
8750  LET  U (7) = (GX-GY) /3 
3760  RETURN 


8160-8240  -  Extract  material 

properties  from  string 
arrays  VAL  converts 
string  to  floating 
point  number 


8540-8580  -  modulus  invariants 


3590-8620  -  Definitions  of  quadratic 
strength  parameters. 

Note  reuse  of  variables 
EX,  ET,  etc. 

8630-8680  -  Convert  from  stress 
space  parameter  to 
strain  space 


8690-8750  -  Failure  parameter 
invariants 
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Figure  1.  Optimization  Trajectory  for  [0/90]  Lai 


42 


rT 


CALLED  FROM  MAIN  UTILITY  ROUTINES 


Figure  3.  Flow  Chart 
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